Computing Counterion Densities at Intermediate Coupling 
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By decomposing the Coulomb interaction into a long distance component appropriate for mean- 
field theory, and a nonmean-field short distance component, we compute the counterion density near 
a charged surface for all values of the counterion coupling parameter. A modified strong-coupling 
expansion that is manifestly finite at all coupling strengths is used to treat the short distance 
component. We find a nonperturbative correction related to the lateral counterion correlations that 
modifies the density at intermediate coupling. 
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The rise of biological physics has rekindled the long- 
standing interest in aqueous electrostatics Q]. Poisson- 
Boltzmann mean-field theory fails to describe a number 
of striking phenomena, such as charge inversion LaLl and 
counterion- mediated attraction 0, la, IS 0> HI El llOj . that 
occur when strong correlations develop between multi- 
valent countcrions. Although there has been some suc- 
cess understanding counterion correlations using both a 
phenomenological Wigner crystal theory |jg, and sys- 
tematic weak-coupling (WC) and strong-coupling 
(SC) 0, 0, 0| expansions, a complete quantitative 
theory spanning the entire range of counterion behav- 
ior is still lacking. This letter introduces a method to 
compute the counterion density at intermediate coupling 
by decomposing the Coulomb interaction into long and 
short distance components in the spirit of the Weeks- 
Chandler- Andersen theory of simple fluids . This de- 
composition not only gives good quantitative agreement 
with simulations, it also provides a natural framework to 
understand both the success of SC expansion, as well as 
the role that lateral correlations play in the counterion 
density. 

Here, I use mean-field theory for the long distance 
interaction, and introduce a modified SC expansion at 
short distances. The traditional SC expansion IS 
problematic because it is formally a virial expansion, and 
one naively expects it to be invalid precisely when the 
counterions are strongly interacting. Nonetheless, nu- 
merical simulations have demonstrated that it not only 
correctly predicts the average counterion density in the 
strong-coupling limit, but also computes the form of the 
corrections 14] . In contrast, the modified SC expansion 
introduced here is manifestly finite in the limit of infinite 
counterion coupling, but recovers the SC corrections at 
large, finite coupling. 

In addition, this decomposition correctly reproduces 
the counterion distribution around a charged surface at 
both strong- and weak- coupling, and agrees very well 
with simulations at intermediate coupling. A test-charge 
theory (TCT) which also computes approximate counte- 
rion densities at intermediate coupling and explains the 



exponential form in the strong coupling limit fails to elu- 
cidate the physics behind the corrections to that limit 
in a clear and satisfactory manner |l6j . Furthermore, I 
find a nonperturbative correction to the density related 
to the lateral counterion correlations that becomes im- 
portant at intermediate coupling. Unlike in SC theory, 
I can unambiguously compute an expression for the free 
energy and show that these lateral correlations play a 
role. 

Consider the primitive model for a charged surface neu- 
tralized by pointlike counterions of the opposite charge 
in a dielectric medium with dielectric constant, e. To 
proceed, introduce a length scale, I, and define V s (r) = 
l B Q 2 e- r ' l /r and Vi(r) = l B Q 2 (l - e-' /£ )A, where l B = 
e 2 / (eksT) is the Bjerrum length and Q is the counterion 
valence. The length I is currently arbitrary and will be 
chosen later to optimize the calculation. The Hamilto- 
nian for ions of charge Qe centered at positions R Q in- 
teracting with a surface of charge density n/(r) = a5(z) 
is given by 

H = \ J d 3 rd 3 r'p(r)[V s (r -r') + Vi{r- r')]p(r'), (1) 

where p(r) = nf(r) /Q — ^ Q 6(r — R Q ). It is understood 
in this expression that ion self-interactions, which arise 
only from V s (r), are to be neglected. The long range 
interaction can be decoupled by introducing a continuous 
field <f> through a Hubbard-Stratonovich transformation, 
resulting in the action S' = S s + Si, where 

-£ [d 3 rn f (v)V s (\r-K a \)/Q (2) 
+ J2V s (\K a -K \), 

£ B — AttIbQ 2 , and 

(3) 
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The counterion positions, R a , are restricted to be over 
the volume of space that can be occupied by the counte- 
rions. 

In the Grand Canonical ensemble for the counterions, 
the partition function is 
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JV^W d 3 R a e- s , (4) 



where the length ft -1 is defined by k 2 /£b = e^/a 3 , a 
is the counterion radius, and p is the chemical poten- 
tial [ill ]. Now define the partial partition function with 
integrals only over the counterion positions, 

oo . „ iV 

z > = E m / II [d 3 R a po(R a )] 
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where p = (k 2 /£ B )e F+l <t> , and F(r) = J d 3 r' V a (\r' 
r\)rtf(r')/Q. This leaves Z oc jT>(f)e~ s , where 
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-lnZ s [p (r)] 
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The average counterion density can be formally com- 
puted as (p(r)) = 5\nZ s /5F{v). 

Since this formulation is exact, the partition function 
is independent of the choice of I. To proceed, make the 
following approximations: (1) the mean- field approxima- 
tion for the long-distance interaction (saddle point in 
</>), and (2) expand the effective potential lnZ s [po] us- 
ing a modified SC-like cluster expansion described be- 
low. Making these approximations, the theory will lose 
its independence on the choice of £, and there will be a 
"best" I whose value gives the closest agreement with the 
full theory. In principle, its value should be determined 
by optimizing the error between a loop expansion on the 
long-distance interaction and perturbative corrections to 
the short-distance expansion. On physical grounds, how- 
ever, we argue that I ~ IbQ 2 , or, in other words, £ is the 
distance that fixed counterions interact with an energy 
of UbT. Counterions at separations larger than this in- 
teract weakly, and mean-field theory is likely valid above 
this length scale. In addition, the distinction between 
short and long length scales should not depend on the 
geometry of the fixed charge distribution, and can there- 
fore only depend on the Bjerrum length and counterion 
valence. For concreteness, I will choose £ — IbQ 2 ■ 

The mean-field approximation, given by 6S/5(j> = 0, 
results in the equation 



V 2 - £ 2 V 4 <£ + £ B (p(r)) = n f (r)£ B /Q, 



(7) 



where I have performed the Wick rotation 4> — > i<p in the 
complex plane for convenience. The SC expansion can be 



reproduced by expanding lnZ s in powers of po(z). How- 
ever, for counterions in the presence of a charged surface 
with charge density cr, each term of this expansion di- 
verges as the coupling constant T — 2ttIbctQ 3 — > oo indi- 
cating that the counterion interactions are not small; this 
divergence can be absorbed by shifting k 2 and utilizing 
overall charge neutrality |l2T ]. 

Instead, expand in powers of dpo(z) = po(z) ~ 
2/(£ B X)5(z), where A = (2tt(tI b Q)- 1 is the Gouy- 
Chapman length for a surface of charge density a. This 
has the property that J dz Spo(z) — due to overall 
charge neutrality, and yields 



Z s = — / d 2 r a dz a Y\_Sp {r a , z a ) 

J ^ 1 ™ 



(8) 
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V s (r a - vp;z a - zp) 



where r a indicates the position of a counterion projected 
to the surface, z a its distance from the surface, and 
Sp {r a , z a ) = (exp[-J2i V s (r a - Ti;z a )])p5po(r a ,Za)- 
Here, (• • -) p is the average taken with respect to the par- 
tition function 



n exp 



i<3 



Vs(ri-r i; 0) 
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The expression (exp[— V s (r a — r^; z a )]) p represents the 
interaction of a charge at coordinates (r Q , z a ) with a layer 
of counterions at positions (rj,Zj = 0). In deriving Eq. 
©, I have assumed (]J a exp[- £\ ^( r « ~ r*;«a)])p ~ 
n„(exp[- £\ Vs{r a ~ ri;z a )]) P , which is true as long as 
the counterions at z > are far enough apart compared 
to L 

Performing a cluster expansion with respect to Spo 
yields lnZ s [<j)} = J d 3 r 5p (r)[l + \ J d 3 r' v 2 (r - 
r')Sp~o(r') + •••], where V2(r) = 1 — exp [— V^(r)]. Terms 
higher than zeroth order vanish in the limit r — > oo as 
the density becomes delta function-like and Spo — ► 0. As 
r — * 0, these corrections also vanish because £ 3 Spo — ► 0, 
and the interactions become predominantly long-ranged. 

It is useful to define po — e^ z ^^ z > with exp[£(r, z)\ = 
(exp[F(z) — J^i ^s(r — ?i] 0)]). The function ((z) can be 
interpreted as the short-distance interaction potential of 
a charge at height z with the charged surface and with a 
layer of counterions at z = 0. Therefore, it encodes the 
response of the z — layer to the presence of a charge 
at some z > 0, and is reminiscent of thcTCT ^ij. One 
difference between £(z) and the TCT, however, is that 
C(z) also depends, at least in principle, on the short- 
range structure of the counterions induced by the short- 
range interaction. To develop a simple approximation 
for £(z) which I will use throughout the remainder of 
the letter, assume that each counterion at z > in- 
teracts with a uniform distribution of charge at z = 
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containing an induced circular correlation hole of radius 
ro = \J £ B X/2 = y/Q/cr. This approximates the size of 
a vacancy in a locally ordered lattice of counterions at 
the surface, which should be valid when T is large. Thus, 

COO - (e/X)[e- z / e -e\ /r o +z2/e }. It is interesting that, for 
z < £, C(z) « T(l - e - r °/ £ ) - z/A is dominated by the 
interaction of the counterions with the bare surface and 
not the z = layer of counterions whose contribution is 
of order — z 2 / (XI) < z/X. 

To lowest order, (p) — po, and the mean- field equation 
for a charged surface now reads 

d 2 z cb - e 2 d 4 z <t> + K 2 e(*)e«»>-*w = e -jfs( z ), (io) 

where cr is the surface charge density. For small z, ({z) — 
F{z) is analytic and can be expanded as a power series 
in (z/£) 2 . On the other hand, F(z) is nonanalytic at 
z = and contributes to the boundary conditions. This 
additional contribution can be disentangled by defining 
<f> = <fi—F. In terms of $, is replaced with C( z ) — -FX Z ) 
and there is an additional source term on the right of Eq. 
(|TU|l of the form —£ B a£ 2 d 2 S(z)/Q. This equation encodes 
two boundary conditions: d z $\%± - £ 2 d 3 z <f>\° ± = £ B a/Q, 
and 9 2 <I>|q+ = £ B v/Q- In terms of the original <f>, the 
boundary conditions are d z <p\ot — 0, and — £ 2 d 3 4>\^_ = 
c£ B /Q- Charge neutrality, J dz k 2 exp[£(z) — <fi(z)] = 
£ B &/Q, is ensured for any solution to Eq. (|10l) . 

Analytical approximations to Eq. I|10[) can be found 
in the WC and SC limits. In the WC limit, I assume the 
solution will decay with characteristic length A. Thus, 
the fourth order derivative is negligible. Since £(r) <c 1, 
cj> has the PB equation form, given by 

(f>(z > 0) = 2 In (l + Kz/Vfj . (11) 

The boundary conditions are satisfied by choosing <fi(z < 
0) = 2(£/X)A(e z / e - 1), where £ 2 A 3 /X 2 - A = 1. This 
requires kX/2 = A. As V — > 0, 4>(z < 0) — > 0, and Eq. 
(|llfl becomes exact. 

In the SC limit, the fourth order term dominates over 
the second order term near the surface. I make the ad- 
ditional assumption that (f> -C 1, and solve 

£ 2 4 Z <P = » n 2 e^~ z /\ (12) 

which has the solution 

(f>{z > 0) = K 2 e c (°) (e~^ A - l) . (13) 

Applying the boundary conditions, <f>(z < 0) = 
-{2/£ 2 )/(l - X 2 /£ 2 )[ex V {z/£) - 1] and k 2 = (2/A 2 )/(l - 
X 2 /£ 2 ) exp[— C(0)]- The exponential SC density is recov- 
ered for large T. This solution also becomes exact as 

r — ► oo. 

For z ^> £, C(z) <C 1 and Eq. (|1(J|) has an approximate 
solution given by Eq. (|11|) . The fourth order derivative 
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FIG. 1: Normalized density difference, n(z) — upb(z), as a 
function of z/X for T — 100 (r = 10 in inset). Solutions to 
Eq. 11011 (solid line) are compared to numerical simulations 
from Ref. (diamonds) and the TCT from Ref. p]| (dashed 
line). 



can be neglected in this limit because k is exponentially 
suppressed by £(0) being large. Therefore, the density at 
large distances is PB-like and is controlled by a renormal- 
ized Gouy-Chapman length, X ren = \/2/n = Aexp[£(0)]. 
This is in agreement with the arguments of Burak et 
al. |lq |. which also exhibits a crossover to a slow decay 
far from the surface. Using ((z), I obtain the estimate 

ln(A ren /A)=r[l- e -VW(^)]. 

Eq. (|10[) has been solved numerically for T = 100 
and r = 10 (inset). Fig. ^plots n(z) — nps(z), where 
n(z) = p(z)£ B X 2 /2 and n PB (z) = 1/(1 + z/X) 2 is the 
normalized Poisson-Boltzmann density. These numerical 
solutions are compared to actual simulation data from 
Ref. 0| (courtesy of A. Moreira) and show quite good 
agreement. Furthermore, Eq. (|10|l outperforms the TCT 
(shown as dashed lines using data provided by Y. Burak 
from Ref. [la|h The nonperturbative function C(z) is 
an important component of this numerical agreement; 
when a virial expansion in po is used to compute In Z s 
(equivalent to setting £(z) = F(z) at lowest order), the 
agreement with the simulation data is only slightly better 
than the TCT and not nearly as good as Fig. ^ A more 
careful evaluation of £(z) is likely to improve these results 
further. To be clear, the value for ro used in Fig. ^ is 
simply an estimate of the real correlation hole size, which 
may differ up to a factor of order one depending on the 
model used. Though there is still very good agreement 
for other reasonable values of ro, ro = y/£ B X/2 seems to 
give the best agreement with simulations. In contrast to 
the r = 10 and 100 results, the density at T = 1 always 
decays faster than the Poisson-Boltzmann density (not 
shown). This is precisely where both the short and long- 
distance expansions in the interpolation scheme attain 
their maximum error, and it is possible that including 
higher order terms may improve this. 

The SC expansion can be reconstructed in this frame- 
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work as an asymptotic expansion around the r = oo limit 
by considering the higher order terms in Spo for lnZ s . 
Here we sketch this result: substituting the conjectured 
asymptotically exact result p = 2e~ z / A / '(^sA 2 ) into (p), 
the first order correction can be computed. It is straight- 
forward to show that this gives only the finite part of the 
first order SC correction as T — * oo, and therefore that 
this correction vanishes in this limit This occurs 

because the delta function in 5po exactly cancels the di- 
vergence from po as T — > oo. Higher order terms will also 
vanish in this limit, and an asymptotic expansion can be 
constructed in powers of 1/r which agrees exactly with 
SC up to first order and which I conjecture also agrees 
at all orders. It is interesting that, although no explicit 
renormalization is necessary in the modified expansion, 
C(z) arises to shift the fugacity k 2 in a similar manner as 
the SC renormalization. A more systematic accounting 
of these corrections is left for future work. 

Despite the quantitative agreement, the modified SC 
expansion suggests a different physical picture of the 
strong-coupling limit: the corrections to the SC limit 
arise from the interactions between only those counte- 
rions that have made excursions away from the wall, as 
measured by Spo- It is clear that the density of these 
excited counterions becomes small for large T, and that 
the SC limit becomes exact. The function ((z) encodes 
the interaction of these excitations with their z = cor- 
relation holes, and becomes important at intermediate 
coupling, once counterions get far enough from the sur- 
face. Interestingly, this correlation hole picture has al- 
ready been described in the context of Wigner crystal-like 
structural correlations for multivalent ions 0, Il8| and in 
the interpretation of the TCT [l^ . 

This scheme also provides a natural framework to com- 
pute the counterion free ene rgy, which is obscured by 
the traditional SC expansion [l2|. This is given by F = 
S (i<p)— pN / (ksT) , where N is the number of counterions, 
cj> is the mean-field long range potential, and the chem- 
ical potential is related to k 2 by p = ksTln^a 3 /£b)- 
Since k 2 depends exponentially on £(()), nonperturbative 
correlations also play a role in the free energy, and subse- 
quently in the interaction between two surfaces at sepa- 
rations where the SC expansion cannot be applied. This 
will be explored in a future publication. 

To summarize, I have computed the counterion den- 
sity around a charged surface using a scheme to decom- 
pose the Coulomb interaction into short and long dis- 
tance components. Each is treated with different ap- 
proximations. For large T, we recover the SC results and 
for small T we recover the WC Poisson-Boltzmann den- 
sity. At intermediate coupling, the model agrees very 



well with the simulation data, and depends crucially on 
a nonperturbative correlation correction whose form we 
have estimated. These correlations also play a role in de- 
termining the renormalized Gouy-Chapman length when 
the densities recovers its Poisson-Boltzman form far from 
the surface. 
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